library(tidyverse)
library(data.table)
library(readxl)
library(GGally)
library(leaflet)
library(raster)
library(mltools)
library(corrplot)
library(highcharter)
Inicialmente, para la lectura de datos, se requieren unos cambios en los nombres de las variables para graficar los mapas posteriormente.
raw_encuesta <- fread("databases/calidad_vida_ok.csv", encoding = "UTF-8") %>%
as_tibble()
# Unificación de las estadísticas de los barrios "ALTAVISTA" y "CABECERA ALTAVISTA"
raw_encuesta <- raw_encuesta %>% mutate(
encuesta_calidad.barrio = case_when(
encuesta_calidad.barrio == "CABECERA ALTAVISTA" ~ "ALTAVISTA",
encuesta_calidad.barrio == "CIUDADELA NUEVO OCCIDENTE" ~ "CABECERA URBANA CORREGIMIENTO SAN CRISTÓBAL",
encuesta_calidad.barrio == "AREA DE EXPANCION SAN CRISTOBAL" ~ "ÁREA DE EXPANSIÓN SAN CRISTÓBAL",
encuesta_calidad.barrio == "CABECERA SAN CRISTÓBAL" ~ "CABECERA URBANA CORREGIMIENTO SAN CRISTÓBAL",
encuesta_calidad.barrio == "SAN CRISTOBAL" ~ "CABECERA URBANA CORREGIMIENTO SAN CRISTÓBAL",
encuesta_calidad.barrio == "PROGRESO Nº 2" ~ "PROGRESO",
TRUE ~ encuesta_calidad.barrio
)
)
# One-hot encoding para variable estudio
variables_categoricas <- raw_encuesta[, c("encuesta_calidad.p_45")] %>% mutate(
encuesta_calidad.p_45 = as.factor(encuesta_calidad.p_45)
) %>% rename(
nivel_estudio = encuesta_calidad.p_45
)
variables_categoricas_encoding <- one_hot(as.data.table(variables_categoricas))
raw_encuesta <- data.frame(raw_encuesta, variables_categoricas_encoding)
La creación de la base de datos se hace con base en las preguntas pertenecientes a la dimensión “Escolaridad”, esta clasificación de las preguntas se hizo a criterio del equipo con la ayuda del informe escrito por la Alcaldía de Medellín durante el planteamiento de la encuesta.
#Agregación de las variables
db_escolaridad <- raw_encuesta %>%
mutate(encuesta_calidad.barrio = str_replace(encuesta_calidad.barrio, "ANDALUCIA", "ANDALUCÍA") %>%
str_replace("Nº 2", "NO.2") %>%
str_replace("Nº 1", "NO.1") %>%
str_replace("Nº 3", "NO.3") %>%
str_replace("AREA EXPANSION", "ÁREA DE EXPANSIÓN") %>%
str_replace("EXPANCION", "EXPANSIÓN") %>%
str_replace("AREA", "ÁREA") %>%
str_replace("BOMBONA", "BOMBONÁ") %>%
str_replace("LA ASOMADERA", "ASOMADERA") %>%
str_replace("BELALCAZAR", "BELALCÁZAR") %>%
str_replace("CALAZANS", "CALASANZ") %>%
str_replace("COLON", "COLÓN") %>%
str_replace("MIRA FLORES", "MIRAFLORES") %>%
str_replace("BARRIO FACULTAD DE MINAS", "FACULTAD DE MINAS") %>%
str_replace("CABECERA SAN ANT DE PR.", "SAN ANTONIO DE PRADO") %>%
str_replace("CARLOS E RESTREPO", "CARLOS E. RESTREPO") %>%
str_replace("URQUITA", "URQUITÁ") %>%
str_replace("LOS CERROS EL VERJEL", "LOS CERROS EL VERGEL") %>%
str_replace("CAYCEDO", "CAICEDO") %>%
str_replace("VALDES", "VALDÉS") %>%
str_replace("CERRO EL VOLADOR", "B. CERRO EL VOLADOR") %>%
str_replace("MOSCU", "MOSCÚ") %>%
str_replace("JOSELA", "JOSÉ LA") %>%
str_replace("JOSE", "JOSÉ") %>%
str_replace("EL YOLOMBO", "YOLOMBO") %>%
str_replace("PIEDRAS BLANCAS", "PIEDRAS BLANCAS - MATASANO") %>%
str_replace("BASILIA", "BRASILIA") %>%
str_replace("VILLA TINA", "VILLATINA") %>%
str_replace("LILIAM", "LILLIAM") %>%
str_replace("BOLIVAR", "BOLÍVAR") %>%
str_replace("CORREGIMIENTO PALMITAS", "PALMITAS SECTOR CENTRAL") %>%
str_replace("INES", "INÉS") %>%
str_replace("FE", "FÉ") %>%
str_replace("LUCIA", "LUCÍA") %>%
str_replace("SABIO", "SAVIO") %>%
str_replace("BERMEJAL- LOS ÁLAMOS", "BERMEJAL-LOS ÁLAMOS") %>%
str_replace("BOLÍVARIANA", "BOLIVARIANA") %>%
str_replace("EL NOGAL - LOS ALMENDROS", "EL NOGAL-LOS ALMENDROS") %>%
str_replace("JUAN XXIII - LA QUIEBRA", "JUAN XXIII LA QUIEBRA") %>%
str_replace("PROGRESO Nº 2", "EL PROGRESO") %>%
str_replace("MARIA", "MARÍA") %>%
str_replace("PLAYÓN", "PLAYON") %>%
str_replace("EL SOCORRO / LA GABRIELA", "EL SOCORRO") %>%
str_replace("FÉRRINI", "FERRINI") %>%
str_replace("LA CANDE LARIA", "LA CANDELARIA") %>%
str_replace("EL PLAYON", "PLAYÓN") %>%
str_replace("IGUANA", "IGUANÁ") %>%
str_replace("MARÍA CANO - CARAMBOLAS", "MARÍA CANO-CARAMBOLAS") %>%
str_replace("DE ABURRA", "DEL ABURRÁ") %>%
str_replace("ALTAVISTA CENTRAL", "ALTAVISTA SECTOR CENTRAL") %>%
str_replace("SECTOR CENTRAL", "CENTRO ADMINISTRATIVO") %>%
str_replace("ALTAVISTA CENTRO ADMINISTRATIVO", "ALTAVISTA SECTOR CENTRAL") %>%
str_replace("SANTA ELENA CENTRO ADMINISTRATIVO", "SANTA ELENA SECTOR CENTRAL") %>%
str_replace("PALMITAS CENTRO ADMINISTRATIVO", "PALMITAS SECTOR CENTRAL") %>%
str_replace("PROGRESO", "EL PROGRESO")
) %>%
group_by(encuesta_calidad.barrio, encuesta_calidad.comuna) %>%
summarize(n = n(), edad_promedio = mean(`encuesta_calidad.p_18`, na.rm = TRUE),
proporcion_lectoescritura = sum(`encuesta_calidad.p_35` == 1, na.rm = TRUE) / sum(!is.na(`encuesta_calidad.p_35`), na.rm = TRUE),
proporcion_escolaridad_actual = sum(`encuesta_calidad.p_36` == 1, na.rm = TRUE) / sum(!is.na(`encuesta_calidad.p_36`), na.rm = TRUE),
estudio_ultimo_ano = sum(`encuesta_calidad.p_37` == 1, na.rm = TRUE) / sum(!is.na(`encuesta_calidad.p_37`), na.rm = TRUE),
proporcion_instituciones_publicas = sum(`encuesta_calidad.p_49` == 1, na.rm = TRUE) / sum(!is.na(`encuesta_calidad.p_49`), na.rm = TRUE),
proporcion_nivel_estudio_0 = sum(`nivel_estudio_0` == 1, na.rm = TRUE) / sum(!is.na(`nivel_estudio_0`), na.rm = TRUE),
proporcion_nivel_estudio_1 = sum(`nivel_estudio_1` == 1, na.rm = TRUE) / sum(!is.na(`nivel_estudio_1`), na.rm = TRUE),
proporcion_nivel_estudio_2 = sum(`nivel_estudio_2` == 1, na.rm = TRUE) / sum(!is.na(`nivel_estudio_2`), na.rm = TRUE),
proporcion_nivel_estudio_3 = sum(`nivel_estudio_3` == 1, na.rm = TRUE) / sum(!is.na(`nivel_estudio_3`), na.rm = TRUE),
proporcion_nivel_estudio_4 = sum(`nivel_estudio_4` == 1, na.rm = TRUE) / sum(!is.na(`nivel_estudio_4`), na.rm = TRUE),
proporcion_nivel_estudio_5 = sum(`nivel_estudio_5` == 1, na.rm = TRUE) / sum(!is.na(`nivel_estudio_5`), na.rm = TRUE),
proporcion_nivel_estudio_6 = sum(`nivel_estudio_6` == 1, na.rm = TRUE) / sum(!is.na(`nivel_estudio_6`), na.rm = TRUE),
proporcion_nivel_estudio_7 = sum(`nivel_estudio_7` == 1, na.rm = TRUE) / sum(!is.na(`nivel_estudio_7`), na.rm = TRUE),
proporcion_nivel_estudio_8 = sum(`nivel_estudio_8` == 1, na.rm = TRUE) / sum(!is.na(`nivel_estudio_8`), na.rm = TRUE),
proporcion_nivel_estudio_9 = sum(`nivel_estudio_9` == 1, na.rm = TRUE) / sum(!is.na(`nivel_estudio_9`), na.rm = TRUE),
proporcion_nivel_estudio_10 = sum(`nivel_estudio_10` == 1, na.rm = TRUE) / sum(!is.na(`nivel_estudio_10`), na.rm = TRUE)
) %>% as_tibble()
# Eliminamos el barrio desconocido
db_escolaridad <- db_escolaridad %>%
filter(encuesta_calidad.barrio != "DESCONOCIDO")
write_excel_csv2(db_escolaridad, "databases/db_escolaridad.csv")
Es importante considerar que la variable “edad_promedio” no necesariamente hace referencia a la dimensión de escolaridad, sin embargo, se quiere realizar un contraste entre la variable “edad_promedio” y las demás variables para extraer información y hacer validaciones empíricas.
summary(db_escolaridad)
## encuesta_calidad.barrio encuesta_calidad.comuna n
## Length:308 Length:308 Min. : 4.0
## Class :character Class :character 1st Qu.: 348.8
## Mode :character Mode :character Median : 862.5
## Mean :1073.2
## 3rd Qu.:1439.0
## Max. :8987.0
## edad_promedio proporcion_lectoescritura proporcion_escolaridad_actual
## Min. :23.00 Min. :0.6667 Min. :0.1000
## 1st Qu.:32.98 1st Qu.:0.8809 1st Qu.:0.2383
## Median :36.20 Median :0.9125 Median :0.2638
## Mean :36.88 Mean :0.9079 Mean :0.2623
## 3rd Qu.:40.69 3rd Qu.:0.9421 3rd Qu.:0.2889
## Max. :58.25 Max. :1.0000 Max. :0.5000
## estudio_ultimo_ano proporcion_instituciones_publicas
## Min. :0.00000 Min. :0.2222
## 1st Qu.:0.01565 1st Qu.:0.7326
## Median :0.02206 Median :0.8948
## Mean :0.02290 Mean :0.8169
## 3rd Qu.:0.02843 3rd Qu.:0.9442
## Max. :0.11111 Max. :1.0000
## proporcion_nivel_estudio_0 proporcion_nivel_estudio_1
## Min. :0.00000 Min. :0.00000
## 1st Qu.:0.07585 1st Qu.:0.06916
## Median :0.12006 Median :0.11016
## Mean :0.12759 Mean :0.11103
## 3rd Qu.:0.16647 3rd Qu.:0.14499
## Max. :0.45455 Max. :0.25397
## proporcion_nivel_estudio_2 proporcion_nivel_estudio_3
## Min. :0.05263 Min. :0.00000
## 1st Qu.:0.19048 1st Qu.:0.06058
## Median :0.25095 Median :0.07895
## Mean :0.23695 Mean :0.07716
## 3rd Qu.:0.28760 3rd Qu.:0.09458
## Max. :0.52381 Max. :0.25000
## proporcion_nivel_estudio_4 proporcion_nivel_estudio_5
## Min. :0.09091 Min. :0.00000
## 1st Qu.:0.20988 1st Qu.:0.03508
## Median :0.24640 Median :0.04964
## Mean :0.24372 Mean :0.04925
## 3rd Qu.:0.27755 3rd Qu.:0.06209
## Max. :0.50000 Max. :0.13333
## proporcion_nivel_estudio_6 proporcion_nivel_estudio_7
## Min. :0.00000 Min. :0.00000
## 1st Qu.:0.01568 1st Qu.:0.02220
## Median :0.03062 Median :0.05376
## Mean :0.03364 Mean :0.09698
## 3rd Qu.:0.04996 3rd Qu.:0.14323
## Max. :0.11111 Max. :0.36842
## proporcion_nivel_estudio_8 proporcion_nivel_estudio_9
## Min. :0.0000000 Min. :0.0000000
## 1st Qu.:0.0004705 1st Qu.:0.0000000
## Median :0.0041754 Median :0.0009129
## Mean :0.0161616 Mean :0.0061769
## 3rd Qu.:0.0198397 3rd Qu.:0.0064606
## Max. :0.1225364 Max. :0.0727763
## proporcion_nivel_estudio_10
## Min. :0.000000
## 1st Qu.:0.000000
## Median :0.000000
## Mean :0.001348
## 3rd Qu.:0.001084
## Max. :0.024691
db_escolaridad %>% select_if(is.numeric) %>%
apply(2, sd)
## n edad_promedio
## 1.019528e+03 5.418267e+00
## proporcion_lectoescritura proporcion_escolaridad_actual
## 4.521456e-02 4.694050e-02
## estudio_ultimo_ano proporcion_instituciones_publicas
## 1.343127e-02 1.811589e-01
## proporcion_nivel_estudio_0 proporcion_nivel_estudio_1
## 6.288825e-02 5.129729e-02
## proporcion_nivel_estudio_2 proporcion_nivel_estudio_3
## 7.227418e-02 2.751268e-02
## proporcion_nivel_estudio_4 proporcion_nivel_estudio_5
## 4.962019e-02 2.082606e-02
## proporcion_nivel_estudio_6 proporcion_nivel_estudio_7
## 2.218821e-02 9.975298e-02
## proporcion_nivel_estudio_8 proporcion_nivel_estudio_9
## 2.558032e-02 1.158639e-02
## proporcion_nivel_estudio_10
## 3.221845e-03
Dada la información sobre las medias arrojada para cada variable, hay que resaltar que la proporción de personas con habilidades lecto-escritoras es alta con un 0.907902, la cual es una cifra positiva pues muestra que en la población no hay analfabetismo crítico.
Por otro lado, las proporción de personas que estudiaron en el último año tambien es baja, pero no se debe concluir que esto es necesariamente negativo. Cabe la posibilidad que varias personas en un determinado barrio sean trabajadores y por tanto no dedican mas tiempo al estudio a partir de cierto momento de su vida.
Ahora centremos nuestra atención en variables específicas del conjunto de datos.
ggpairs(data = db_escolaridad %>% dplyr::select(edad_promedio, proporcion_escolaridad_actual, proporcion_lectoescritura),
mapping = aes(alpha = 0.7))
De la gráfica anterior, se debe resaltar significativamente las correlaciones entre las variables “edad_promedio” con las variables “proporcion_lectoescritura” y “proporcion_escolaridad_actual”. Note que la “edad_promedio” esta correlacionada de manera directa con la variable “proporcion_lectoescritura” y de manera inversa con la variable “proporcion_escolaridad_actual”. Esta información que se obtiene coincide con el conocimiento epírico, pues por un lado, las personas con mayor edad tienen mas probabilidad de tener como mínimo educación primaria, obteniendo así mas habilidades lecto-escritoras, mientras que por otro lado a medida que las personas envejecen dejan de educarse debido a que pasan al entorno laboral.
ggplot(data = db_escolaridad,
mapping = aes(x = proporcion_instituciones_publicas)) +
geom_density()
Por otro lado, analizando la gráfica anterior, se puede notar que en los barrios, las personas tienen tendencia a educarse en establecimientos de índole pública debido a su fácil acceso.
Ahora centraremos nuestro análisis en la proporcion de personas cuyo último nivel de estudio fue primaria, secundaria y universidad. A continuación interpretaremos la relación entre estas proporciones y el tipo de institución pública donde se realizó el último estudio.
ggpairs(data = db_escolaridad %>%
dplyr::select(proporcion_instituciones_publicas,
proporcion_nivel_estudio_7),
mapping = aes(alpha = 0.7))
En la gráfica anterior se puede ver un hallazgo importante, y es que en los barrios donde las personas tienden a estudiar en instituciones públicas, hay una menor tendencia a optar por estudios universitarios, mientras que los barrios que tienen tendencia a realizar estudios en instituciones privadas, tienen mayor proporción de personas que optan tener estudios universitarios. Este fenómeno puede analizarse desde una perspectiva transversal a la condicion económica de los barrios, pues si en un barrio hay más tendencia a educarse en instituciones privadas, significa que hay un poder adquisitivo mucho más alto, lo que permite tener oportunidades económicas para realizar estudios universitarios; mientras que por otro lado, el acceso a la educación pública para estudios escolares es económicamente mas viable, por tanto en los barrios con bajos ingresos la educación primaria y secundaria se puede adquirir con facilidad, pero ingresar a instituciones de educación superior es mucho más retador, pues las universidades privadas no son una opción común debido a las dificultades económicas, y sumado a esto la educación superior pública tiene mecanismos de selección para estudiantes aspirantes, lo cual reduce la posibilidad de ingresar a la educación universitaria en estos barrios.
Ahora si tomamos las proporciones de los niveles de estudio y extraemos la media puede tenerse una visión general de la escolaridad en Medellín.
medias_nivel_estudio <- apply(db_escolaridad %>% dplyr::select( proporcion_nivel_estudio_0,
proporcion_nivel_estudio_1,
proporcion_nivel_estudio_2,
proporcion_nivel_estudio_3,
proporcion_nivel_estudio_4,
proporcion_nivel_estudio_5,
proporcion_nivel_estudio_6,
proporcion_nivel_estudio_7,
proporcion_nivel_estudio_8,
proporcion_nivel_estudio_9,
proporcion_nivel_estudio_10), 2, mean)
ggplot(data = data.frame(nivel_estudio = factor(0:10), medias_nivel_estudio)) +
geom_bar(mapping = aes(x = nivel_estudio, weight = medias_nivel_estudio))
Note entonces que las personas de los barrios tienden a tener estudios de secundaria y media técnica como último estudio realizado, mientras que los estudios universitarios tienen una proporcion media muy baja respecto a los demas niveles de estudio.
Inicialmente escalaremos las variables que se usaran en los algoritmos de agrupamiento.
db_escolaridad_escalado <- db_escolaridad %>%
mutate_at(vars(-n, -encuesta_calidad.barrio, -encuesta_calidad.comuna), scale)
# Selección de las columnas para el agrupamiento
datos_escalados_agrupamiento <- db_escolaridad_escalado %>%
dplyr::select(-n, -encuesta_calidad.barrio, -encuesta_calidad.comuna)
pca <- prcomp(datos_escalados_agrupamiento)
summary(pca)
## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6
## Standard deviation 2.840 1.4951 1.10595 0.96429 0.94599 0.76194
## Proportion of Variance 0.504 0.1397 0.07645 0.05812 0.05593 0.03628
## Cumulative Proportion 0.504 0.6437 0.72012 0.77824 0.83417 0.87045
## PC7 PC8 PC9 PC10 PC11 PC12
## Standard deviation 0.73055 0.63054 0.54745 0.49933 0.45841 0.41797
## Proportion of Variance 0.03336 0.02485 0.01873 0.01558 0.01313 0.01092
## Cumulative Proportion 0.90381 0.92866 0.94739 0.96297 0.97611 0.98703
## PC13 PC14 PC15 PC16
## Standard deviation 0.31929 0.23779 0.22158 2.654e-16
## Proportion of Variance 0.00637 0.00353 0.00307 0.000e+00
## Cumulative Proportion 0.99340 0.99693 1.00000 1.000e+00
Ahora, se mostrará la varianza explicada por los componentes mediante el calculo de la suma acumulada del porcentaje explicado.
porcentaje_exp <- 100 * (pca$sdev ^ 2) / sum(pca$sdev ^ 2)
suma_acumulada_porcentaje <- cumsum(porcentaje_exp)
ggplot(data = data.frame(componente = factor(1:ncol(datos_escalados_agrupamiento)), suma_acumulada_porcentaje),
aes(x = componente, y = suma_acumulada_porcentaje, group = 1)) +
geom_point() +
geom_line() +
geom_hline(yintercept = 90, color = "red") +
geom_vline(xintercept = 7)
En la gráfica anterior, la linea horizontal roja corresponde a una varianza explicada acumulada de un 90%, por tanto se puede ver que al tomar las primeras 7 componentes principales se obtiene una varianza explicada de 90.3809835%, lo cuál se usará para hacer la reducción de dimensionalidad.
Primero proyectaremos los datos respecto a las componentes principales.
datos_proyectados <- t(t(pca$rotation) %*% t(datos_escalados_agrupamiento))
Ahora seleccionaremos las primeras 7 componentes principales que explican mas del 90% de la varianza.
datos_agrupamiento <- datos_proyectados[, 1:7] %>% as_tibble()
El primer método que se usara para hacer el agrupamiento sera K-means. Para esto, realizaremos una gráfica de codo que nos dará indicios del número de grupos que debemos tomar dentro de este algoritmo.
k_maximo <- 15
wss <- sapply(1:k_maximo,
function(k) {
kmeans(datos_agrupamiento, k, nstart = 10, iter.max = 15 )$tot.withinss
})
ggplot(data = data.frame(k = factor(1:k_maximo), wss),
aes(x = k, y = wss, group = 1)) +
geom_point() +
geom_line() +
xlab("Número de grupos [K]") +
ylab("Total de suma de cuadrados intra-grupos")
k_optimo_kmeans <- 5
Note entonces que una selección de \(K\) entre 5 y 7 es razonable, ya que la variación de la suma de cuadrados intra-grupos no ofrece un cambio significativo. Además se considera contraproducente hacer una selección muy grande de grupos en términos de la interpretabilidad del resultado ofrecido por el algoritmo de agrupamiento, pues se considera mas util tener pocos grupos que sean altamente diferenciadores de la ciudad de Medellín contrario a tener una variedad de grupos altamente dispersos donde no se pueda identificar un patrón de agrupamiento claro.
Se realiza la ejecución del algoritmo K-means con la cantidad de grupos extraidas del análisis anterior.
agrupamiento_kmeans <- kmeans(datos_agrupamiento, centers = k_optimo_kmeans, nstart = 10, iter.max = 15 )
Al colorear el mapa de los barrios de Medellín de acuerdo a los grupos se obtiene la siguiente gráfica.
political <- shapefile("Barrio_Vereda/Barrio_Vereda.shp")
Encoding(political@data$NOMBRE) <- "UTF-8"
political$NOMBRE <- political$NOMBRE %>% toupper() %>% str_replace("DE MESA", "DE MESA")
grupos_barrios <- data.frame(barrio_nombre = db_escolaridad$encuesta_calidad.barrio, comuna_nombre = db_escolaridad$encuesta_calidad.comuna, grupo = agrupamiento_kmeans$cluster)
nombres_mapa <- data.frame(nombre_barrio = political$NOMBRE)
vector_nombres = c()
vector_grupos = c()
for(nombre_mapa in nombres_mapa$nombre_barrio) {
grupo <- grupos_barrios[grupos_barrios$barrio_nombre == nombre_mapa, 3][1]
vector_nombres <- c(vector_nombres, nombre_mapa)
vector_grupos <- c(vector_grupos, grupo)
}
mapa_grupos <- data.frame(
nombre_barrio = vector_nombres,
grupo = vector_grupos
)
factpal <- colorFactor(rainbow(k_optimo_kmeans), mapa_grupos$grupo)
leaflet(data = political) %>%
addTiles() %>%
addProviderTiles(providers$OpenStreetMap) %>%
addPolygons(fill = TRUE, stroke = TRUE, weight = 2, color = ~factpal(mapa_grupos$grupo),
label = as.character(paste(political$NOMBRE, " - Grupo ", mapa_grupos$grupo)),
popup = as.character(paste(political$NOMBRE, " - Grupo ", mapa_grupos$grupo))) %>%
addLegend("bottomright", colors = "#03F", labels = "Barrios y veredas")